vquantile <- function(X,W=rep(1/NROW(X),NROW(X))){
  # This function solves the semi discrete optimal transport problem
  # and obtains the vector quantile using otm2D, a function from
  # the Geogram library based on the Multi-scale approach of Merigot [2011].
  #
  # The vector quantile obtained here is a "power diagram": a list of vertices of 
  # convex polygons (cells) that symbolize the mapping from [0,1]x[0,1] to the data.
  # Output of this function is intended to pass to "ILF.R". 
  #
  # Output: pwd        = power diagram representing the vector quantile from OT
  #         weight     = The weights of the power diagram to recreate pwd
  #         constraint = max difference between areas of the cells and constraint
  #                      to roughly test if constraints are satisfied
  #         flag       = counts the number of empty cells
  #                      if any cells return empty, it is a fail of convergence, 
  #                      repeat algorithm.
  #         time       = time elapsed
  #
  # Inputs: X = bivariate allocation (data)
  #         W = sampling weights of X, default is equal weights (1/NROW(X))
  #
  #Summary: Step 1- Ensure there are no duplicated vectors, otherwise algorithm fails.
  #         Step 2- Normalize mean of X to one.
  #         Step 3- Calculate the vector quantile. If empty cells occur, repeat.
  #
  start <- Sys.time()
  df <- data.frame(X,W)
  colnames(df) <- c("X1","X2","weight") 
  # Step 1
  df <- aggregate(weight ~ X1 + X2, df, sum) # aggregate based on duplicates
  # Step 2
  W <- df$weight/sum(df$weight)
  mu <- cbind(df[,1]%*%W, df[,2]%*%W)
  X <- cbind(df[,1]/mu[1],df[,2]/mu[2]) # weighted mean normalized to 1
  # Step 3
  Y <- scaling.min.max(X) # scale data to make compatible with OT problem
  redo <- 1
  s <- 0
  while (redo == 1){
    if (s > 9){return("Issue: Empty cells, check data.")}
    s <- s + 1
    redo <- 0
    tryCatch({
      ot <- otm2D(Y,W) # Calculate weights of power diagram
      pwd <- power_diagram(Y[,1],Y[,2],ot$weights) # Calculate power diagram
      # Collect diagnostic information on if the constraint is satisfied
      # and if there are no empty cells
      constraint <- sapply(which(is.na(pwd$cells)== FALSE), function(j){polyarea(pwd$cells[[j]][,1], pwd$cells[[j]][,2])})
      constraint <- max(sqrt(sum((constraint - W)^2)))
      empty <- length(which(is.na(pwd$cells)))
      if ((empty > 0) || (constraint > 10^(-4))){redo <- 1}
    },error = function(e){redo <<- 1})
  }
  end <- Sys.time()
  return(list(pwd = pwd,weight=ot$weights,constraint = constraint,flags = empty, time = end - start))
}